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Detailed protocol 


SUGAR-seq protocol for tumor infiltrating lymphocytes (TILs) 


FACS-sorted TILs were incubated with a mixture of biotinylated and non-biotinylated L-Pha (1:5 ratio), to prevent signal saturation, at a concentration of 1 g/ml 
for 20 minutes on ice. After three washes with cold phosphate-buffered saline (PBS), cells were stained with cell hashing antibodies and CITE-seq (cellular 
indexing of transcriptomes and epitopes by sequencing) antibodies as previously described (7). An anti-biotin oligo-conjugated CITE-seq antibody was used 
to convert the L-Pha signal into a sequenceable readout. “Stained” and washed cells were counted, brought to ~1000 cells/ul, and loaded onto the 10x 
chromium instrument (10x Genomics, Pleasanton, CA, USA) to generate single-cell gel beads in emulsion and capture/barcode cells. For 5’ multiomics 
including gene expression, TCR-sequencing and hashing/CITE-seq, we used the 10x Single Cell V(D)J Kit with Feature Barcoding (enabled for 5' gene 
expression, TCR, and feature barcoding for cell surface protein) protocol following the manufacturer’s instructions. TCR libraries were prepared using the 10x 
Chromium Single Cell V(D)J Enrichment Kit, mouse T cell. 5’ hashtag oligos (HTO)/ADT (one library) was prepared using the Chromium Single Cell 5' Feature 
Barcode Library Kit and indexed using the Chromium i7 Multiplex Kit N, Set A. Cell Ranger reads were aligned to the mm10 reference genome; cellular 
barcodes were demultiplexed, and unique molecular identifiers and antibody capture (ADT and HTO) were quantified using 10x Genomics’ Cell Ranger 
software (version 3.1.0). Cell barcodes containing RNA or antibody counts from cells from more than one sample (intersample doublets) were identified using 
Seurat’s HTODemux function. Barcodes containing counts from more than one cell within the same sample (intrasample doublets) were identified using the 
Scrublet Python package (version 0.2.1) (2). A cutoff of more than one median absolute deviation value above the median Scrublet score was chosen. Cells 
identified as either type of doublet were removed from the analysis. Relative TF activity in each cell was estimated using the single-cell regulatory network 
inference and clustering method (3). Area-under-curve (AUC) scores for TFs were calculated on the basis of inferred gene regulatory networks using the 
pySCENIC Python package (version 0.9.19). Gene expression and antibody count matrices were processed in R (version 3.6.1) using the Seurat R package 
(version 3.1.0) (4).RNA transcript counts for barcodes identified as cells by Cell Ranger were normalized using sctransform via Seurat’s SC Transform function. 
The counts were first transformed with no covariates in the sctransform model; then, cell cycle phase scores were estimated using Seurat’s CellC ycleScoring 
function with mouse homologs of the cell cycle gene sets provided by Seurat. The sctransform normalization was then rerun with the cell cycle phase scores and 
the percentage of raw RNA counts belonging to mitochondrial genes for each cell as variables to be regressed out in the model. ADT counts were log1p 
transformed and scaled (zscore) or normalized using centered log ratio transformation with either method, producing a similar result. Principal components 
analysis was then performed on the sctransform-scaled RNA expression values for genes with residual variance in the sctransform model greater than 1.3.A 
shared nearest neighbor (SNN) network was calculated using the top 10 principal components using the FindNeighbors function with k-nearest neighbors set to 
50 and cosine distance metric. The SNN network was then used to identify cell populations using the FindClusters function using the Louvain algorithm with 
resolution parameter 0.6. UMAP values were also calculated using the RunUMAP function with the top 10 principal components as input and parameters 
n.neighbors = 50 and metric = “cosine.” Monocle3 was used for pseudotime analysis using default setting. 
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